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Abstract 

An exactly solvable model for the rewiring dynamics of weighted, directed networks 
is introduced. Simulations indicate that the model exhibits two types of condensation: 
(i) a phase in which, for each node, a finite fraction of its total out-strength condenses 
onto a single link; (ii) a phase in which a finite fraction of the total weight in the 
system is directed into a single node. A virtue of the model is that its dynamics can be 
mapped onto those of a zero-range process with many species of interacting particles 
— an exactly solvable model of particles hopping between the sites of a lattice. This 
mapping, which is described in detail, guides the analysis of the steady state of the 
network model and leads to theoretical predictions for the conditions under which the 
different types of condensation may be observed. A further advantage of the mapping 
is that, by exploiting what is known about exactly solvable generalisations of the zero- 
range process, one can infer a number of generalisations of the network model and 
dynamics which remain exactly solvable. 

1 Introduction 

Networks, and in particular weighted networks, have a long history in the social, natural and 
engineering sciences. The many and varied examples of networks range from systems such 
as transportation networks to ecological, biochemical and social networks |I1E|- Basically, 
a network is a graph containing nodes connected by links (or edges) — the network is said 
to be 'weighted' when weights are assigned to the links. The interest in networks within the 
physics community was initially with regard to the statistical properties of the connectivity 
structure of unweighted networks (i.e., networks where links are restricted to have weight 
one or zero only, indicating a connection or no connection respectively). Studies focussed, for 
example, on the degree distribution, i.e., the distribution of the number of links connected 
to a node, and the clustering coefficient which probes how trios of nodes are correlated. 
Subsequently simple models for growing networks exhibiting interesting statistical properties 
were proposed [3]. 



Recently, weighted networks have also become of interest to this community. These net- 
works may be defined through their adjacency matrix Wij, which gives the weight of the 
link from node i to node j. The links may be directed or undirected corresponding to an 
asymmetric or symmetric adjacency matrix. Data are available for many real weighted net- 
works, of either directed or undirected kind, leading to studies of the statistical properties of 
Scientific Collaboration Networks (SCN) and the World-wide Airport Network (WAN) 
ini Hj • Typical quantities used to describe the topology of general weighted networks are: 
the distribution of in-strength and out-strength, which are defined as the total weight going 
into or out of a node due to all the connected links; weighted clustering coefficients; gener- 
alisations of the idea of minimal spanning trees to weighted networks |H1 IH] • An interesting 
aspect of weighted networks is how the weights of the links are correlated to the topology 
of the network: the weights may be determined solely by the global topology, they may be 
independent, or they may be correlated in some more complicated way. For example, in 
Transportation Networks, where the weight represents the fiow through a link, the weights 
are determined solely by the topology; in the SCN it appears that the degree of a node and 
weights of its links are uncorrelated; in the WAN, correlations between the node degree and 
the weights of its links lead the strength to grow faster than the degree. 

Attention has recently focussed on constructing simple models of 'growing' weighted 
networks, with the aim being to investigate the resulting degree and weight distributions and 
the relation between the two. Initial work had the weight of links determined by the degree 
of the node they were attached to UHl UH 112] • To loosen this coupling between weight 
and degree, a model with dynamical evolution of weights during growth was introduced in 
IT^ ITK] . Moreover, models of growing networks in which the dynamics depend on the 
weights of the links rather than the degree of the nodes were proposed in For a review 
of recent work on weighted networks see [7] . 

An alternative class of evolving network models, which is the focus of this work, is that 
of 'rewiring' or 'equilibrium' networks In these models, usually considered for 

unweighted networks, the number of nodes is fixed and the dynamics involves the rewiring of 
links between nodes, although the effect of creating and destroying links has been examined 
recently jTHj. For rewiring dynamics, one is interested in the steady-state properties of an 
ensemble of networks. Under certain conditions on the rewiring rules [201 one can obtain a 
condensed phase where one node (the condensate) attracts a finite fraction of the available 
links. Thus there is a single extensively connected node while the other nodes exhibit a power- 
law background in the degree distribution. Hence the network is scale free if one discounts 
the condensate. The transition is reminiscent of the 'gelation' transition in growing networks 
from the scale free phase to a phase with a dominant 'hub' pil I22[ 1^. 

These rewiring networks may be related to interacting particle systems such as the zero- 
range process [21] for which condensation transitions have been widely studied — for a recent 
review see [201 • In the zero-range process, particles hop between the sites of a lattice with a 
rate that depends on the number of particles at the site of departure. The idea behind the 
mapping is to identify the nodes of the network with the sites of the lattice model and the 
links of the network become the particles — the rewiring of links corresponds to particles 
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hopping from site to site. At the simplest level, this mapping is approximate in that the 
number of particles at a site represents only the degree of the corresponding node, but does 
not determine the other nodes to which a node is attached. However, the steady states of 
the two systems have the same form. 

The aim of the present work is to introduce a model for the rewiring dynamics of a 
directed, weighted network. By virtue of an exact mapping between the network model 
and a set of coupled interacting particle systems, we are able to analyse exactly conditions 
on the rewiring rates under which the network will exist in some kind of condensed phase. 
Specifically, we define the network model in Section |21 the key features of which are: 

1. The network is directed. 

2. The out-strength of each node is conserved under the dynamics but the in-strength 
evolves. 

3. The dynamics are governed by the weights of links and the in-strength of the nodes. 

In addition, the weights are integer variables, although this constraint may be relaxed (as 
we discuss in Sectional). Thus we define a weighted network model with a fixed number of 
nodes where strong nodes and links with large weight may tend to attract more weight. 

In Section in| we show how this model enjoys a mapping to an interacting particle system 
known as a multi-species zero-range process [SEIEZI- The idea is that, since the out-strength 
of each node is conserved, the dynamics at each node is basically the exchange of weight 
between the links coming out of that node. Thus each node may be related to an interacting 
particle system in which the exchange of particles between sites represents the exchange 
of weight between links. However, since the dynamics at each node also depends on the 
in-strength, these interacting particle systems are coupled. This mapping of the weighted 
network to a set of coupled interacting particle systems contrasts with the previous mappings 
of an unweighted network to a single interacting particle system discussed in [T71 ITHll^ . 

The results of simulations, presented in Section El indicate two distinct types of con- 
densation within this model: the first is when, for any given node, a finite fraction of the 
out-strength condenses onto one weight out of that node; the second is when a finite fraction 
of the total weight in the system is directed onto one particular node. The first scenario 
gives a realisation of a transition in the 'disparity', the idea of 'disparity', discussed in the 
context of internet traffic |2H1 and simple growing weighted network models i2£j , being that 
one weight from a node is dominant over the others from that node. The second scenario 
is related to the usual condensation transition occurring in unweighted rewiring networks, 
wherein a single 'hub' acquires a finite fraction of the total number of links in the system. 

We exploit the mapping to the interacting particle system in Section |3J in order to solve 
and analyse steady state properties of the network model exactly. In particular, we derive 
conditions on the rewiring rates which lead to the two types of condensation observed in 
simulations. The mapping also enables one to identify a number of possible generalisations 
of the network model and its rewiring dynamics, which we discuss in Section Finally, we 
conclude in Section IHl 
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2 Model 



We study a model of a dynamically evolving, directed network of L nodes with weighted 
links. Since the weights are integers we will denote the weight associated with the link from 
node k into node / by an integer riki > (instead of the usual Wki)- When Uki = 0, this 
represents no link from k to I. We depict such a network with L = 4 nodes in Figure ^ 

The dynamics of the network are such that one unit of the weight associated with the 
link pointing from node k into node / is rewired to point into another randomly selected 
node I' with a rate Uk{ni), where Ui = nu, . . . ,nLi- In general, this rate is a function of the 
weights associated with all of the links pointing into node I, and it depends on the source 
node k: we will specialise to a rate of the form 

L 

Ukiiii) = u%'^ki)u''C^nki) , (1) 
fc=i 

which depends on the weight of the link being rewired, through the function w'^lriki), and the 
total in-strength of node /, through the function u'^{Ylk=i ^ki)- So for the example in Figure 
^ the rate at which the one unit of weight for the link connecting node C to node A is rewired 
to another node is mcI'^aa, "'^ba, ''^CA, ""-da) = m''(3)-u'^(10). After the rewiring, ricA = 2 and 
any of ucb, ncc or ncn have increased by one. These dynamics conserve the out-strength 
= Ylf=i ^ki which is the total weight of all links pointing out of node k; the total weight 
of all links pointing into node /, i.e., the in-strength Xi = J2k=i^ki, is not conserved. For 
simplicity we will take Mk = M for all k. Later, in Section we will consider several 
generalisations of the model and dynamics. In particular we consider rewiring dynamics 
which depend on the node to which a link is being rewired. This allows one to make a 
connection with a form of preferential attachment in which links are preferentially rewired 
to nodes with large in-strengths ^H]- We proceed with the dynamics defined above in order 
to make a clear connection to well-studied interacting particle systems. We stress that either 
alternative for the dynamics yields the same behaviour and transitions. 

2.1 Mapping to a system of interacting particles 

To make the connection with interacting particle systems, we use the adjacency matrix. A, 
to encode the network: the value of element (fc, /) of the matrix is n^/, the weight of the link 
pointing out of node k into node /, as shown for example in FigureHJ Note that here, because 
the network is directed, A ^ A^ . The idea is to represent the adjacency matrix by a lattice, 
where site (/c, /) of the lattice represents element (fc, /) of the adjacency matrix. The value 
of the element (fc,/), Uku is then understood as the number of particles at site (/c, /) of the 
lattice. The lattice is therefore two-dimensional, since it represents the elements of a matrix, 
even though the network model may typically be defined in a fully connected geometry (i.e. 
a link pointing into one node is rewired to point into any other randomly selected node). 

The network dynamics then causes the elements of the adjacency matrix to change. 
Thus the corresponding particle configurations of the lattice model also evolve according 
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Figure 1: A picture of a weighted directed network with L = 4 nodes, labelled A, B, C and 
D. The direction and weight of each link is shown. The corresponding adjacency matrix is 
also shown. 



to dynamical rules, rules which have their origin in the dynamics of the network model. 
For the network dynamics described above the corresponding interacting particle system is 
a generalisation of the much studied zero-range process ,20,, for which the steady state is 
exactly solvable. 

2.2 Model as a ZRP with L species of particles 

The interacting particle system is defined on a square lattice (since the adjacency matrix is 
a square matrix) with L x L sites. On this lattice Uki labels the number of particles at the 
site located at row k = 1, . . . , L, column I = 1, . . . , L. There are M particles in each row, a 
total of LM in the system. Under the network dynamics described in Section |2l both M and 
L are conserved. Rewiring single units of the weight of a link connecting node k to node / 
corresponds, in the interacting particle system, to a particle in row k and column / hopping 
to another randomly selected site within the same row with a rate Uk{n^), given by (P). 

We remark that these are the dynamics of a ZRP with L species of particles, labelled 
k = 1, . . . ,L. Each site / = 1, . . . , L of a one- dimensional lattice contains n^i particles 
of species k. A particle of species k then hops with a rate Uk{nj) to any other site on 
the lattice, i.e. a rate which depends upon the number of particles of each species at the 
departure site I. Thus we identify rows of the two-dimensional lattice with different particle 
species (so the total number of particles of each species is conserved) and the columns of the 
two-dimensional lattice are identified with sites of a one-dimensional lattice, as illustrated 
in Figure El This mapping is useful because much is known about the exact steady states 
of zero-range processes and many of its generalisations. The model defined above can be 
solved exactly in the steady state: the analysis which follows exploits this knowledge and 
one expects that many variations of the network model will also be exactly solvable, those 
which are can be inferred from the corresponding zero-range process. 

Although the model we consider can be thought of as a ZRP with L species of particles, 
we will proceed with the terminology of the rows and columns of a two-dimensional lattice. 
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Figure 2: A picture of the interacting particle systems which represent the network depicted 
in Figure n On the left-hand side is the two-dimensional lattice model with the occupation 
numbers Uki entered in row k and column / of the lattice. On the right-hand side is the cor- 
responding multi-species ZRP (here there are four species) — note that there is no ordering 
amongst the particles at a particular site. 



since this retains closer contact with the adjacency matrix of the original network model. 



3 Simulations 

In this section we present results from simulations of the model that show the two types of 
condensation: site condensation, where a finite fraction of the particles in each row condenses 
onto a single randomly located site of the row; and column condensation, where a finite 
fraction of all the particles in the system condenses into a single column. In the network 
context this corresponds to a finite fraction of the weight from each node being contained in 
a single link from that node; and a finite fraction of the weight pointing from all nodes to a 
single node, respectively. We also present simulations which exhibit characteristics of both 
types of condensation which are interesting from a networks perspective. 

The model is simulated using a simple Monte Carlo algorithm. At each time step the 
following update procedure is followed: 

1. A site of the lattice {k, I) is chosen at random. 

2. If the site is occupied then a particle will be removed from this site with probability 
u{ri])St, where u is the hop rate ([Q) and 6t is the time interval, chosen such that the 
probability of a hop occurring is always less than or equal to one. 

3. If a particle was removed, then a second site {k, m) (within the same row) is chosen at 
random and a particle is added to this site. 

All simulations were run on lattices of size 100 x 100, i.e. on networks with 100 nodes, 
and for 0(10'') Monte Carlo Sweeps. The first half of each run was used to relax to the 
steady state and the second half to measure probability distributions. Distributions were 
measured for 
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• the number of particles at a site (link weight) 

• the number of particles in a column (node in-strength) 

• the number of occupied sites in a column and row (in- and out-degree respectively). 

Typical configurations of the lattice were also output. 

In order to compare with theory (see Section E}, the following hop rates were chosen to 
show the two condensation types. Recall that we consider a hop rate for particles from site 
{k, I) of the form (^. For site condensation we make the choice 

u%n) = (1 + ¥/n) and u%n) = 1 , (2) 

with b"^ = A and for n > 0. For column condensation we make the choice 

u^in) = < I ~^ / — ' g^^^ u^in) = 1 , (3) 

^ ^ l + b^L/n for n>L, \ j y \ j 

with 5^ = 1.05. 

3.1 Site Condensation 

The simulations for the site condensation case were run with 175 particles in each row 
starting from a random initial configuration. The site condensation can be seen clearly from 
the site distribution (black circles, Figure 01 (a)). The peak at around n = 150 represents 
the site condensates, it has an area of order 1/L and there are sites indicating that L 
sites have occupations of around 150 particles. As the number of particles in each row is 
fixed at 175, there can be at most one such site in a row and an occupation of this size 
represents an appreciable fraction of the total number of particles available to a site, hence 
implying a condensation in the thermodynamic limit. Thus we have a single condensed site 
in each row. Note also that the site distribution has an apparently power-law background 
before the condensate peak, this is reminiscent of the behaviour of a single-species ZRP [20] . 
indeed as the column attraction has effectively been switched off, the system is a collection 
of single-species ZRPs. The column distribution (grey crosses. Figure El (a)) shows that, in 
the absence of any coupling between the rows, the condensed sites are randomly distributed 
among the columns — as one would expect. The first peak in the distribution at around 
n = 50 shows that some columns do not contain any condensed sites, just those from the 
power-law background. The second, third, . . . peaks correspond to columns with one, two, 
. . . condensed sites in them, plus a background from the rest of the sites. This is exactly 
what one would expect if the condensed sites were randomly distributed in the columns and 
a random sampling of the site distribution to create a column distribution agrees closely 
with the simulation. The typical configuration (Figure |3)) bears this out, condensates can be 
seen on individual sites, but no column ordering can be discerned. Although not shown here, 
systems with weak column attraction showed similar behaviour, i.e. as though no column 
attraction was present. 
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The degree distributions for the site condensation case (Figure El (b)) appear to be bino- 
mial in nature: a binomial distribution constructed from the measured probability of a site 
having zero occupation matches very closely the data shown. Thus in the network context, 
for the site condensation we have connectivity similar to that of a random graph (which 
also has a binomial connectivity distribution). In addition to this we have a condensed link 
weight from each site with a power-law background, and an in-strength which is a random 
sum of link weights. 

3.2 Column Condensation 

The simulations for the column condensation case were run with 1000 particles in each row. 
The initial condition was taken with all particles on the same site. The column condensation 
can be seen clearly from the column distribution (grey crosses. Figure (a)). The peak at 
around n = 85000 corresponds to the condensate, it has an area of order 1/L and there are L 
columns, indicating that a single column has an occupation of around 85000, an appreciable 
fraction of the total number of particles in the system Lx M = 100 x 1000. The other, lower 
peak at around n = 100 represents the columns which do not contain a condensate. The site 
distribution (black circles. Figure (b)) also shows condensation onto sites in the peak at 
around n = 850. If the column interaction were to be switched off leaving L uncoupled single- 
species ZRPs, it is known that the chosen u^{n) (= constant) would not give condensation 
in this fully connected homogeneous system, see for example [20] • Thus in this case it is 
the column interaction that has induced a site condensation. Note that both the column 
and site distributions have apparent power-law pieces, with equal or close exponents. The 
theory presented in Section E] allows one to construct the critical column distribution at the 
transition point — shown as the dotted line in Figure El (a). This fits the non-condensate 
background part of the measured column distribution very well. The typical configuration 
of the system. Figure IHl shows the column condensate comprises site condensates all in the 
same column. 

The out-degree distribution for the column condensation case (grey triangles, FigureEKb)) 
is similar to the site condensation case, being distributed according to a binomial distribution. 
However, the in-degree distribution (Figure El (b), black squares) is significantly different, 
in particular it has a broader tail and also a single site which is connected to all others, 
as shown by the lone data point at degree L with probability 1/L {L = 100 in this case). 
Thus in the network context the out-connectivity resembles that of a random graph, but the 
in-degree distribution displays a broader tail and there is a single node to which all others 
point — the node corresponding to the condensed column. Furthermore, all links pointing 
to this node hold a finite fraction of the weight available and the node holds a finite fraction 
of the total in-strength available to it. Broadly-tailed degree distributions are often observed 
in real networks, weighted and unweighted alike. While the tail observed in these simula- 
tions is not as broad as many observed, it is at least broader than that of the equivalent 
random graph. Also apart from the condensed pieces the in-strength and weight distribu- 
tions display power-law tails, again something that has been observed in many real networks. 
Thus, as with equilibrium networks [21 Ej , certain distributions may take a power-law form 



8 



(a) 



10" 



10 



1" 



site (weight) distribution 
column (in-strength) distribution 



10" 



10- 



10" 



10"' 



10 




10 

10"' 



10 



10 



10 
n 



10 



10 



(b) 



10" 



10 



10' 



fi -4 



10" 



10"^^ 



10" 
10"' 



10" 



a 
a 

H 



6 
6 



6 
A 



□ in-degree distribution 
A out-degree distribution 



10 



10 



10^ 



n 



Figure 3: Distributions from simulation of the model with the hop-rate chosen to display site 
condensation, (a) Site and column occupation distributions, corresponding to the distribu- 
tions of weight among links and in-strength among nodes in the network context respectively, 
(b) In- and out-degree distributions. 
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Figure 4: Typical configuration of tfie system, output from simulation of the model with the 
hop-rate chosen to display site condensation. 



at a critical point. Power-law strength and weight distributions and non-power-law degree 
distributions have been observed for a traffic network in 

As mentioned at the beginning of this subsection, we employed a fully condensed initial 
condition. The reason is that, for a random initial configuration, the dynamics slows down 
as the system reaches a state in which several columns contain a large number of particles; 
the single column condensate is attained only on prohibitively long timescales. 

3.3 Further behaviour 

In the preceding sections the somewhat extreme cases of the absence of site or column 
attraction in the presence of the other (m* or u'^ set equal to one) were considered. This was 
to allow for a direct comparison with theory. Cases where some forms of non-trivial attraction 
are present in both and are difficult to compare directly with theory. However, such 
systems do display interesting behaviour and so in this section we present some numerical 
data from simulations of such a system. 

Simulations were run on a system with 175 particles in each row and the hop-rate given 
by: u'^{n) = 1 + h'^ / L for n < L, u'^{n) = 1 + b'^/n for n > L and u'^{n) = 1 + b'^/n, with 
W = 16 and 6* = 2.5. If were set equal to one, then no condensation would occur at 
this particle number |2Q|. If u"^ were set equal to one, then no condensation would take 
place at any finite particle number (see Section 0]). Results from simulations show that 
this system, with both site and column attraction present in forms that separately show no 
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Figure 5: Distributions from simulation of the model with the hop-rate chosen to display 
column condensation, (a) Site and column occupation distributions, corresponding to the 
distribution of weight among links and in-strength among nodes in the network context 
respectively, (b) In- and out-degree distributions. 
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Figure 6: Typical configuration of tfie system output from simulation of the model with the 
hop-rate chosen to display column condensation. 



condensation behaviour, does show something that could be interpreted as a condensation- 
like phenomenon. The following results were obtained from a random initial configuration. 

In Figure m (a), it can be seen that the site occupation distribution (black circles) has 
a decaying part and a peak at around n = 150 and is generally reminiscent of a condensed 
system. The column distribution (Figure d (a) grey crosses) also shows a large n peak, at 
around n = 3500, although it is much broader than one usually associated with a clear 
condensation. The typical configuration data in Figure |H1 indicate that the high n bump of 
the column distribution represents more than one highly occupied column, each composed 
of many site condensates. Without a direct comparison with theory available, it is difficult 
to interpret this broad peak: it could be a finite-size effect that will not be seen in larger 
systems, but on the other hand it could be a true condensation that is being adversely 
affected by the finite size of the system. In either case it is interesting from a network point 
of view as real- world networks are often of finite size. 

The degree distributions of this system are also interesting. They both have a binomial 
form for low degree, but the in-degree distribution departs from this at high degree — it has 
a small secondary peak implying that several sites are highly connected, but a completely 
connected site to which all others point, as arises in the column condensation case, is absent. 
Thus in the network context we have a network with connectivity much like that of a random 
graph, except for the existence of several well-connected 'hub' nodes which would not be 
present in the random case. These hub nodes also have many high weight links pointing to 
them, giving them a large in-strength. 
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For an initial configuration with all particles located on a single site, we obtain distribu- 
tions of the same form as above and configurations typically contain more than one highly 
occupied column. 



4 Theory 

In this section we present the exact steady state of the model and exploit this solution in 
order to understand theoretically the condensation transitions observed in simulations. 



4.1 Steady state 

First, we give the steady state, the proof of which is given in an appendix. The steady state 
probabilities -P({z^z}), for finding the system in a configuration {ui}, are given by the simple 
form 



(4) 



1=1 



i.e. a product over columns. Here, the functions /(nj depend on a single column (i.e. the 
in-links of node /), and for the rates (0), they are given by 



(5) 



k=l 



where we are using Xi = Ylk=i ^ki to represent the total number of particles in a column 
(the node in-strength), and where 

X 

nx)=l[u%z)-\ (6) 

1=1 

for a = s,c, where /°(0) = 1. The normalisation Z^^m in © is given by 



" = E n 

{"■ki} '=1 



k=l 



l[6{J2nki-M) 



(7) 



k=l 1=1 



The 5-function in this normalisation enforces particle conservation within each row (which 
is the conservation of node out-strength in the network): it is this 5-function that induces 
correlations between different columns in the steady state. 



4.2 Condensation theory 

Now, we exploit the exact steady state in order to understand theoretically the condensation 
transitions observed in simulations. Ideally, one would wish to demonstrate condensation in 
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Figure 7: Distributions from simulation of the model with the hop-rate chosen such that 
without the site or column attraction elements no condensation would be present, (a) Site 
and column occupation distributions, corresponding to the distributions of weight among 
hnks and in-strength among nodes in the network context respectively, (b) In- and out- 
degree distributions. 
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Figure 8: Typical configuration of the system output from simulation of the model with the 
hop-rate chosen such that without the site or column attraction elements no condensation 
would be present. 



the site- and column-distributions of particles, p'^{n) and p^{n). The expressions for these 
distributions involve the normalisation Zl^m given by ((Zj), which can be thought of as a 
canonical partition function since M is fixed. However, as we now outline, it turns out to 
be simplest to work within the grand canonical ensemble in which the particle number is 
allowed to fiuctuate. 



4.2.1 Grand canonical ensemble 



We introduce the grand canonical partition function in the usual way, by defining fugacities 
{zk} so that we can replace the canonical partition function ((Zj) by 



oo L L 
k 

1=1 



{nfe,=0} k=l 



k=l 



where the fugacities are chosen to ensure that, on average, each row contains the proper 
number of particles M^: 

Zk^ = 2_^nki = Mk. (9) 

Here, the overbar indicates an average taken in the grand canonical ensemble. Since we are 
interested in the case where Mk = M for all values of /c, the Ihs of this equation must be 
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independent of k therefore Zk = z for all values of k; in this case, the fugacity z is chosen to 
determine the total number of particles in the system. After a little rearrangement, (jS)) can 
be written 



n 

1=1 



L \ L 

\k=l nj^i / k=l 

-1 L 



k=l 



(10) 

(11) 



where in going from the first to the second line we use the fact that the I subscript on the n^i 
plays no role — each column makes the same contribution to — and the square bracket 
can simply be raised to the power of L. Now, X = ra^. By taking derivatives of Zl with 
respect to z, one finds that the fugacity must be chosen to satisfy 



M 



dlnF{z) 
dz 



where we have defined the function 



F(z) = 5]r(x)nr(^.) 



yrik 



(12) 



(13) 



k=l 



The equation (fT^ is the key result of this subsection. Condensation typically arises 
when one is unable to satisfy p2j) above some critical density; under these circumstances, 
the grand canonical ensemble is no longer valid and the change of the partition function from 
the grand canonical form below the critical density to some other form above it signals the 
phase transition. We now discuss conditions under which (fT^ can be satisfied in order to 
illustrate how one can theoretically understand and predict the condensed phases observed 
in simulations. 



4.2.2 Condensation 

It is straightforward to show that the function F{z) is a smoothly increasing function of its 
argument z. Let us assume that we have chosen hop rates such that the sums which determine 
F{z) have a radius of convergence a. The convergence properties of F{z) and its derivative 
determine whether or not the system undergoes condensation. These convergence properties 
are determined by the form of the function f{n) which in turn provides the conditions on the 
rewiring rates under which the system condenses. We now present two very simple examples 
to illustrate the two types of condensation. 

First, we consider conditions to observe site condensation. We take /^(a;) = 1, which is 
the case when u'^[x) = 1 for all x, hence 



F{z) 



(14) 



n=0 
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Hence, from ()12|). the fugacity is determined by 



where p = M/L is the density of particles in each row. When both the numerator and 
denominator approach finite values as z approaches the radius of convergence a there exists 
a finite critical density pc above which (fT^ cannot be satisfied. This signals condensation: for 
p > Pc, each row contains (p — pc)L 'excess' particles which condense in a single, randomly 
located site within each row. In order that pc is finite, one requires that /^(n) decays 
asymptotically like /^(n) ~ for large n where b > 2. This in turn implies that the 

hop rates M*(n) must decay asymptotically as 

u'{n) a{l + b/n) , (16) 

again, where 6 > 2, if one is to observe site condensation in each row at finite density. In 
the network model, this implies that the rate of rewiring a node must decay more slowly 
than 2/n, where n is the weight of the link pointing out of the node being rewired. These 
considerations guide the choice of hop rates © discussed in Section ITTl 

Now, we consider conditions to observe column condensation. To deal with the sum over 
particle configurations in F{z), it turns out to be convenient to sum over the variable X 
and introduce a 5-function to ensure that X represents the number of particles in a column. 
Thus F[z) may be written in the form 

oo L L 

{„j.}X=0 k=l k=l 

where the k subscripts refer to the rows within a single column. To proceed, we take 
f^in) = 1, which is the case when u^{n) = 1 for all n. With this choice, the sum over {rik} 
is straightforward to perform — the 5-function constraint supplies a combinatorial factor 
associated with the number of ways of adding L integers to get X — and one finds 

X ^ ' 

In this case, the rhs of ()12|) represents the average number of particles in a column, X. The 
asymptotics of f^iX) determine the convergence properties of F{z) (fTH|) . In order to deduce 
the required asymptotics, we expand the binomial factor for X ^ L, using the approximation 

L-l+X\ X^-i , , 

(19) 



L-1 J (L-1)!' 

for large X. To observe column condensation, the rhs of |T2|l must converge to a value 
Xc = 0{L) as 2; — i> a. This is the case when /'^(X) decays asymptotically like /'^(n) ~ 
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a~^X~'*' with A > 1 + L. When these conditions are satisfied, the system condenses above a 
critical density pc = Xc/ L in such a way that every site in a single, randomly located column 
typically contains (X — X^) particles. A hop rate which yields f^iX) with the required 
asymptotics is 

u'{X) ~ a(l + cL/X) , (20) 

where c > 1, which motivates the choice of hop rates discussed in Section IX^ 

Thus these two simple examples, /'^(x) = 1 and /''(n) = 1, illustrate site and column 
condensation respectively. We would like to stress that in both site and column condensation, 
while only the asymptotics determine whether or not the system condenses, the actual value 
of the critical density depends on the details of the form of the rewiring rates for all values 
of their arguments, not just the asymptotics. 

We now indicate how one can use the grand canonical ensemble to predict the background 
particle distribution in the condensed phase. We consider the case f ''^{n) = 1, for which 
the probability that a column contains exactly X particles, is given by 

f{X) = ^ . (21) 

This equation holds throughout the low density phase. In the condensed phase, the grand 
canonical ensemble breaks down, however (PT|) with z = a correctly reproduces the form of 
the background distribution of column occupation numbers. For the rates Q, ()2H) cannot 
be written in a convenient form at 2; = 1, however it can easily be evaluated numerically 
for a given value of b'^. The result of such a computation, with b'^ = 1.05, is compared with 
simulation in Figure El (a). 



5 Generalisations 

In this section we consider some generalisations of the network model which, in the steady 
state, are still given by a factorised form. The aim is to illustrate that the class of network 
models with factorised steady states extends to a wide variety of rewiring dynamics and so, 
by making the connection between network dynamics and interacting particle systems, we 
are rewarded with a versatile approach to the analysis of network properties. 



5.1 General dependence of the rewiring rates on departure and 
destination nodes 

We begin by considering three generalisations of the rewiring rates: 

1. General dependence of the rewiring rates on the weights of links rather than the spe- 
cialised form ((T)). 

2. Rewiring rates which depend on the weights of the links pointing into both the node 
from which the link is being removed, and the weights of the links pointing into the 
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node into which the hnk is being rewired — this allows one to consider preferential 
attachment. 

3. Heterogeneity in the rewiring dynamics such that one can consider rates which differ 
depending on the source, departure and destination nodes of the link. 

With these generalisations, one unit of the weight associated with the link pointing from 
node k into node / is rewired to point into another randomly selected node /' with a rate 
''J'kiilh) iki'ijh')- The heterogeneity described in the above point 3 enters through the /c, /, and 
/' superscripts on u and t. The corresponding interacting particle system has been discussed 
in the literature in the context of 'Urn' models (discussed in ^31j) and 'Misanthrope' processes 
(discussed in [20]) in a fully connected geometry, in which particle hop rates depend on the 
numbers of particles at both the site of departure and the destination site, but now defined 
with L species of particles. 

In the steady state, this model can still be expressed in the factorised form 

L 

P{{rh}) = Z-L^MX{fiilk). (22) 
1=1 

(though note fi{n^ now carries a sub-script /) provided the hop rates satisfy the constraint 
Ukijlh] nk'i-'^) tkijui; Uki-l) _ Uk'iilh'^ nki-l) tkniHi; Uk'i-l) 

for all pairs of rows k and k' . We have introduced the notation n^;nki—l which represents 
the configuration but with the term n^i replaced by riki—l. When the hop rates satisfy 
this constraint, the functions fi{ni) can be written 

L 

Mm) = n 

k=l 

having chosen /(O, ... ,0) = 1. This can be written in a number of different forms due to 
fl23|) — the symmetry in fi{ni) is obscured within this constraint. A proof of this steady 
state is presented in Appendix B. 

5.1.1 Preferential attachment 

To illustrate how one might exploit this steady state, we consider a network model with 
preferential attachment rewiring dynamics. Thus we choose Ukiini) = 1 for all k and /, and 

L 

tki'im') = t'inu')t^C^nu') , (25) 

k=l 

which defines the rate a link pointing from node k into node / is rewired to point into node 
This rate therefore is a function of the weight of the link pointing from the source node k 



Y\ ^kijO, ■ ■ . ,0,i-l,nk+ii, . . .,nLi) 
UkiiO, ■ ■ ■ ,0,i,nk+ii, ■ ■ ■ ,nu) 



(24) 
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into the destination node /' before the rewiring event, and a function of the total in-strength 
of the destination node. It satisfies the constraint (j23|) and one finds that the functions f{ni) 
are given by (0), where 

X 

nx) = l[t%z-l), (26) 

i=l 

for a = s,c, where /"(O) = 1. In order that this model will exhibit the behaviour observed 
in simulations, one can follow the analysis of Section 14.2.21 to deduce that if the hop rates 
decay asymptotically like 

f(n)~l-6/n, (27) 

with b > 2, then one observes site condensation above some finite critical system density, 
and if 



f (X) ~ 1 - X/X , (28) 

for X ^ L, with A > L + 1, then one observes row condensation above some finite critical 
system density. 



5.2 Further generalisations 

There are a number of further ways one can choose to redefine the model and still obtain 
a factorised steady state. One way is to consider geometries other than the fully connected 
one, such as rewiring dynamics in which the weight being rewired from node I is always 
rewired to a 'neighbouring' node I + 1. All the models considered so far also have factorised 
steady states in this geometry although there may exist extra constraints in certain cases. 
Generalisations to more complicated geometries are also possible, and have been discussed 
in the context of the ZRP in 

We can also relax the property of integer weights and consider the case where the weights 
are continuous variables and an arbitrary amount of the weight is allowed to be rewired. The 
approach is similar to the way one generalises the zero-range process to continuous masses 

It is straightforward to allow each node to have a different total out-strength, but more 
complicated generalisations whereby these out-strengths are not conserved may also be pos- 
sible. Non-conservation in the context of the ZRP has been discussed in jl9t l2Uj. Another 
straightforward modification of the model is to prohibit links which point into the same node 
they point out of, i.e., restrict to Ukk = 0. A further possibility is to allow negative weights, 
as well as positive, which corresponds to a generalisation of the so-called 'Bricklayers' model' 
[5^ — positive and negative weights have been considered in a rewiring social network in 

A point about the factorised steady states we would like to emphasise is that one is free 
to choose any form for /(n^) and then infer the rewiring rates from a recursion such as ()B.3jl . 
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6 Conclusion 



In summary we have introduced a dynamical model of a weighted, directed network wherein 
the out-strength of each node is conserved but the in-strength is not. The rewiring dynamics 
of this model can be mapped onto a multi-species zero-range process for which the steady 
state is exactly solvable. 

Through numerical simulations and theoretical analysis we have identified two conden- 
sation transitions in the steady state. In the network context they correspond to a disparity 
transition where for each node a single link contains a finite fraction of the out-strength, and 
the more familiar condensation transition where the in-strength of a single node captures a 
finite fraction of the out-strengths for all nodes in the system. These transitions demonstrate 
some of the varied behaviour that is possible in weighted directed networks. 

Within the multi-species ZRP picture these transitions correspond to, in the first case, 
condensation of all species at independent sites or, in the second case, a collective condensa- 
tion of all species onto the same site. From our knowledge of the multi-species ZRP we are 
able to identify a number of generalisations of the model and dynamics which preserve the 
exactly solvable steady state. For example these correspond in the network context to contin- 
uous weights, preferential attachment dynamics and node fitness. It would be interesting to 
explore further the different possible behaviours afforded by these and other generalisations. 
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Appendices 

Appendix A. Proof of steady state (j4l) to (EI) 

The steady state (jHIHl) can be demonstrated by noting that, since the network model is defined 
in a fully connected geometry (i.e. a link can be rewired from a node to any other node in 
the network), the steady state probabilities must satisfy detailed balance with respect to the 
rewiring dynamics. (In Sectional we comment on how one generalises to other geometries.) 
The detailed balance condition implies that 

u'{nM+lX{Xi+l)P{{ni}; Uki+l, n^^-l) = u'{nk„^)u%X^)P{{n^}) , (A.l) 

where we have introduced the notation {nj; n^^+l, rijtm— 1 to represent a particle configu- 
ration with Uki replaced by Uki+l and Ukm replaced by Ukm—^- In this balance equation, 
the Ihs represents a particle hop from site kl to site km and the rhs represents the reverse 
process. We now regard a solution of the form ()4l5p as an ansatz, substitution of which into 
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the balance equation ()A.1|) yields 

after cancelling common factors and rearranging slightly. This is satisfied for 

= rix) , (A.3) 
for a = s, c, which is iterated to give (jHl). This proves the steady state (jH) to (jHl). 

Appendix B. Proof of steady state (|22D to (l24h 



The steady state (j22j) to ()24|) can be demonstrated by asking that the steady state proba- 
bilities satisfy detailed balance with respect to the rewiring dynamics, hence 

for each k = 1, . . . ,L and all pairs / 7^ Inserting the steady state ((221), cancelling common 
factors and rearranging yields 

Ukiim'^riM+l) fi{ni;nki+l) _ Ukvijii,) fvijh') 2) 

tkiirn) fiilk) tki\lh']nkv-l) fv{n^,]nkv-l) ' 

again, for each k = 1, . . . , L. Both sides of this equation must be equal to a constant which 
we choose, without loss of generality, to be equal to one, whereby we obtain the recursion 

fiilk) = r~\ .^'(^i' ^ki-l) , (B.3) 



which can be iterated to obtain ()24|). This recursion also implies the constraint on the choice 
of the hop rates: one can obtain two expressions for fi{ni) in terms of fi{ni;nki—l,nk'i—^)', 
one comes by applying the recursion ()B.3|) first to k, then to k'; the other comes by applying 
the recursion to k' before k. Performing these operations yields the constraint (j23j) . 

References 

[1] R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 47, 74 (2002) 

[2] S. N. Dorogovtsev and J. F. F. Mendes 2003 Evolution of Networks (OUP, Oxford) 

[3] A.-L. Barabasi and R. Albert, Science 286, 509 (1999) 

[4] M. E. J. Newman, Proc. Natl. Acad. Sci. USA 98, 404 (2001) 

[5] M. E. J. Newman, Phys. Rev. E 70, 056131 (2004) 



22 



[6] R. Guimera and L. A. N. Amaral, Eur. Phys. J. 38, 381 (2004) 

[7] M. Barthelemy, A. Barrat, R. Pastor-Satorras and A. Vespignani, Physica A 346, 34 
(2005) 

J. D. Noh and H. Rieger, Phys. Rev. E 66, 066127 (2002) 



[9] P. J. Macdonald, E. Almaas, A.-L. Barabasi, |cond-m at/0405688| 

[10] S. H. Yook, H. Jeong, A.-L. Barabasi and Y. Tu, Phys. Rev. Lett. 86, 5835 (2001) 

[11] T. Antal and P. L. Krapivsky, Phys. Rev. E 71, 026103 (2005) 

[12] E. Almaas, P. L. Krapivsky and S. Redner, Phys. Rev. E 71, 036124 (2005) 

[13] A. Barrat, M. Barthelemy and A. Vespignani, Phys. Rev. Lett. 92, 228701 (2004) 

[14] A. Barrat, M. Barthelemy and A. Vespignani, Phys. Rev. E 70, 066149 (2004) 

[15] A. Barrat, M. Barthelemy and A. Vespignani, Lect. Notes Comput. Sci. 56, 3243 (2004) 

[16] S. N. Dorogotsev and J. F. F. Mendes, |cond-mat/0408343| 



[17] Z. Burda, J. D. Correia and A. Krzywicki, Phys. Rev. E 64, 046118 (2001) 

[18] S. N. Dorogovtsev, J. F. F. Mendes and A. N. Samukhin, Nucl. Phys. B 666, 396 (2003) 



[19] A. G. Angel, M. R. Evans, E. Levine and D. Mukamel, |cond-mat/0503487 



[20] M. R. Evans and T. Hanney, J. Phys. A 38, R195 (2005) 

[21] P. L. Krapivsky, S. Redner and F. Leyvraz, Phys. Rev. Lett. 85, 4629 (2000) 

[22] P. L. Krapivsky and S. Redner, Phys. Rev. E 63, 066123 (2001) 

[23] G. Bianconi and A.-L. Barabasi, Phys. Rev. Lett. 86, 5632 (2001) 

[24] M. R. Evans, Braz. J. Phys. 30, 42 (2000) 

[25] B. Hu, X.-Y. Jiang, J.-F. Ding, Y.-B. Xie and B.-H. Wang, |HOTld^mat/0408125 

[26] M. R. Evans and T. Hanney, J. Phys. A 36, L441 (2003) 

[27] T. Hanney and M. R. Evans, Phys. Rev. E 69, 016107 (2004) 

[28] M. Barthelemy, B. Gondran and E. Guichard, Physica A 319, 633 (2003) 

[29] G. Bianconi, |cond-mat/04i2'399| 



[30] A. De Montis, M. Barthelemy, A. Chessa and A. Vespignani, cond-ma t/05071Q6| 

23 



[31] C. Godreche and J. M. Luck, J. Phys.:Condens. Matter 14, 1601 (2002) 

[32] M. R. Evans, S. N. Majumdar and R. K. P. Zia J. Pliys. A: Math. Gen. 37, L275 (2004) 
[33] M. Balazs, Annales de I'lnstitut Henri Poincare - PR 39, 639 (2003) 



24 



